Euler
ExamplesWe create an adjacency matrix for Penney's game. For the game, see
It contains the probabilities to go from a triplet of coin tosses (x,y,z) to another one (y,z,w) with an additional coin toss. The triplets are coded in dual form, i.e., (0,0,0) is stored as triplet 0, (0,0,1) as triplet 1 and so on. Note, the EMT starts its vectors with index 1, however.
In the matrix, M[i,j] contains the probability to go from triplet j to triplet i in one coin toss. Most of the transitions are impossible.
>function makeM () ...
M = zeros(8,8);
for i1=0 to 1;
for i2=0 to 1;
for i3=0 to 1;
for j1=0 to 1;
for j2=0 to 1;
for j3=0 to 1;
i=4*i1+2*i2+i3;
j=4*j1+2*j2+j3;
if i2==j1 and i3==j2 then
M[j+1,i+1]=1/2;
endif;
end;
end;
end;
end;
end;
end;
return M;
endfunction
>M=makeM; >shortest M
0.5 0 0 0 0.5 0 0 0
0.5 0 0 0 0.5 0 0 0
0 0.5 0 0 0 0.5 0 0
0 0.5 0 0 0 0.5 0 0
0 0 0.5 0 0 0 0.5 0
0 0 0.5 0 0 0 0.5 0
0 0 0 0.5 0 0 0 0.5
0 0 0 0.5 0 0 0 0.5
Let us simulate a random walk, starting at triplet (0,0,0).
>p=zeros(8)'; p[1]=1; p'
[1, 0, 0, 0, 0, 0, 0, 0]
>p=M.p; p'
[0.5, 0.5, 0, 0, 0, 0, 0, 0]
>p=M.p; p'
[0.25, 0.25, 0.25, 0.25, 0, 0, 0, 0]
>p=M.p; p'
[0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125]
After only three steps we have equal chances to reach every other triplet.
Now we walk through the triplets randomly, but stop at two specific triplets T and S. We want to learn how many walks stop at T and how many walks stop at S.
>M=makeM;
We need to cut off all connections from S and T. In the example, we take
S = (0,0,1) i.e. triplet in column 2 T = (0,1,0) i.e. triplet in column 3
>M[,2]=0; M[2,2]=1; // S is (0,0,1) >M[,3]=0; M[3,3]=1; // T is (0,1,0) >shortest M
0.5 0 0 0 0.5 0 0 0
0.5 1 0 0 0.5 0 0 0
0 0 1 0 0 0.5 0 0
0 0 0 0 0 0.5 0 0
0 0 0 0 0 0 0.5 0
0 0 0 0 0 0 0.5 0
0 0 0 0.5 0 0 0 0.5
0 0 0 0.5 0 0 0 0.5
We start with equal probability at all triplets.
>p=ones(8)'/8; p'
[0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125]
>p=M.p; p'
[0.125, 0.25, 0.1875, 0.0625, 0.0625, 0.0625, 0.125, 0.125]
>p=M.p; p'
[0.09375, 0.34375, 0.21875, 0.03125, 0.0625, 0.0625, 0.09375, 0.09375]
Now we do 100 steps.
>loop 1 to 100; p=M.p; end; >p'
[0, 0.666667, 0.333333, 0, 0, 0, 0, 0]
This shows that triplet (0,0,1) is superior ot triplet (0,1,0) with a winning probability of 2:1.
This time,
S = (1,0,0) i.e. triplet in column 5 T = (0,0,0) i.e. triplet in column 1
>M51=makeM; >M51[,5]=0; M51[5,5]=1; >M51[,1]=0; M51[1,1]=1; >shortest M51
1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0.5 0 0 0 0.5 0 0
0 0.5 0 0 0 0.5 0 0
0 0 0.5 0 1 0 0.5 0
0 0 0.5 0 0 0 0.5 0
0 0 0 0.5 0 0 0 0.5
0 0 0 0.5 0 0 0 0.5
>p=ones(8)'/8; p'
[0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125]
>loop 1 to 1000; p=M51.p; end; >p'
[0.125, 0, 0, 0, 0.875, 0, 0, 0]
>p[5]/p[1]
7
We learn that (1,0,0) wins against (0,0,0) with a chance of 7:1.
Let us denote the probability that a triplet S is reached before a triplet T when we start at knot i by w(i). Then we have the properties
- w(S) = 1 - w(T) = 0 - w(i) is the average w(j) for all j that can be reached from i.
Thus
I.e., w is an eigenvalue of Mij'.
>V=eigenspace(M51',1)
1 0
0 0.377964
0 0.377964
0 0.377964
0 0.377964
0 0.377964
0 0.377964
0 0.377964
We want an eigenvector which is 0 at the T=5 and 1 at S=1.
>v = V[,2]/V[5,2]
0
1
1
1
1
1
1
1
Then we start on one triplet randomly and get for the average win.
>sum(v'/8)
0.875
We can make that automatic and compute all probabilities that S appears before T if started at a random triplet.
>function getw (M,i,j) ... if i==j then return 0; endif; Mt=M; Mt[,i]=0; Mt[i,i]=1; Mt[,j]=0; Mt[j,j]=1; V = eigenspace(Mt',1); c = V[[i,j]]\[1;0]; return sum((V.c)')/cols(Mt); endfunction
>M = makeM();
Let us try our simulated examples.
>fraction getw(M,2,3), fraction getw(M,5,1)
2/3 7/8
Now we can compute a matrix W. The entry at position (i,j) is the probability to get to triplet number i before triplet number j if we start at a random triplet.
>function solvePG () ...
M = makeM();
W=zeros(8,8);
for i=1 to 8;
for j=1 to 8;
W[i,j]=getw(M,i,j);
end;
end;
return W;
endfunction
>W=solvePG(); >fracformat(7); W, defformat;
0 1/2 2/5 2/5 1/8 5/12 3/10 1/2 1/2 0 2/3 2/3 1/4 5/8 1/2 7/10 3/5 1/3 0 1/2 1/2 1/2 3/8 7/12 3/5 1/3 1/2 0 1/2 1/2 3/4 7/8 7/8 3/4 1/2 1/2 0 1/2 1/3 3/5 7/12 3/8 1/2 1/2 1/2 0 1/3 3/5 7/10 1/2 5/8 1/4 2/3 2/3 0 1/2 1/2 3/10 5/12 1/8 2/5 2/5 1/2 0
The result is that in each column there is one number bigger than 1/2. If player A selects triplet number j, we select triplet number i such that W(i,j)>1/2.
We want to print the better triplet is human readable form.
>function printtriplet (i) ... s = ")"; i=i-1; if mod(i,2) then s=",1"+s; else s=",0"+s; endif; i = floor(i/2); if mod(i,2) then s=",1"+s; else s=",0"+s; endif; i = floor(i/2); if mod(i,2) then s="(1"+s; else s="(0"+s; endif; return s; endfunction
>function printPG (W) ... for j=1 to 8; w=W[,j]'; e=extrema(w); printtriplet(j) + " -> " + printtriplet(e[4]), end; endfunction
We get the well known solution of the game.
>printPG(W)
(0,0,0) -> (1,0,0) (0,0,1) -> (1,0,0) (0,1,0) -> (0,0,1) (0,1,1) -> (0,0,1) (1,0,0) -> (1,1,0) (1,0,1) -> (1,1,0) (1,1,0) -> (0,1,1) (1,1,1) -> (0,1,1)